
clear
capture log close
set more off


******* cpi 2
. clear

. use "bls_prices_2.dta"

* normalize cpi to equal one in base years

foreach num in 2005  {
	gen norm`num' = cpi if year==`num'
	sort year
	replace norm`num' = norm`num'[_n-1] if norm`num'==.
	gsort -year
	replace norm`num' = norm`num'[_n-1] if norm`num'==.
	gen cpi_`num' = cpi/norm`num'
	
}
* save
keep year cpi_*
sort year
save cpi_allyears_2005, replace


clear
use cps_00003
sort year serial pernum
merge m:1 year using cpi_allyears_2005
*keep if _merge==3
drop _merge

* calculate real income

gen hhincome_2005 = hhincome/cpi_2005

** keep household head only
keep if relate==101
	
* Census division
gen division = 1 if region==11
replace division = 2 if region==12
replace division = 3 if region==21
replace division = 4 if region==22
replace division = 5 if region==31
replace division = 6 if region==32
replace division = 7 if region==33
replace division = 8 if region==41
replace division = 9 if region==42

save income_b2005.dta,replace

clear
use income_b2005.dta

forvalues i=1990/2020 {
use income_b2005.dta

* create real income groups
	gen hhfaminc_2005 = 1 if hhincome_2005<2500
	replace hhfaminc_2005 = 2 if hhincome_2005>=2500 & hhincome_2005<5000
	replace hhfaminc_2005 = 3 if hhincome_2005>=5000 & hhincome_2005<7500
	
	replace hhfaminc_2005 = 4 if hhincome_2005>=7500 & hhincome_2005<10000
	replace hhfaminc_2005 = 5 if hhincome_2005>=10000 & hhincome_2005<15000
	replace hhfaminc_2005 = 6 if hhincome_2005>=15000 & hhincome_2005<20000
	replace hhfaminc_2005 = 7 if hhincome_2005>=20000 & hhincome_2005<25000
	replace hhfaminc_2005 = 8 if hhincome_2005>=25000 & hhincome_2005<30000
	
	replace hhfaminc_2005 = 9 if hhincome_2005>=30000 & hhincome_2005<35000
	replace hhfaminc_2005 = 10 if hhincome_2005>=35000 & hhincome_2005<40000
	
	replace hhfaminc_2005 = 11 if hhincome_2005>=40000 & hhincome_2005<45000
	replace hhfaminc_2005 = 12 if hhincome_2005>=45000 & hhincome_2005<50000
	replace hhfaminc_2005 = 13 if hhincome_2005>=50000 & hhincome_2005<55000
	replace hhfaminc_2005 = 14 if hhincome_2005>=55000 & hhincome_2005<60000
	replace hhfaminc_2005 = 15 if hhincome_2005>=60000 & hhincome_2005<65000
	replace hhfaminc_2005 = 16 if hhincome_2005>=65000 & hhincome_2005<70000
	replace hhfaminc_2005 = 17 if hhincome_2005>=70000 & hhincome_2005<75000
	replace hhfaminc_2005 = 18 if hhincome_2005>=75000 & hhincome_2005<80000
	replace hhfaminc_2005 = 19 if hhincome_2005>=80000 & hhincome_2005<85000
	replace hhfaminc_2005 = 20 if hhincome_2005>=85000 & hhincome_2005<90000
	replace hhfaminc_2005 = 21 if hhincome_2005>=90000 & hhincome_2005<95000
	
	replace hhfaminc_2005 = 22 if hhincome_2005>=95000 & hhincome_2005<100000
	replace hhfaminc_2005 = 23 if hhincome_2005>=100000 & hhincome_2005<120000
	replace hhfaminc_2005 = 24 if hhincome_2005>=120000


keep if year==`i'

collapse (sum) asecfwt, by(hhfaminc_2005)
rename asecfwt hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save `i'_cps_2005,replace
}






*** predictefd kwh

clear
forval i=1990/2020  {
clear
use recs_2005
capture drop if kwh==999999

bysort moneypy: egen sumhh=total(nweight)
bysort moneypy: gen weight=nweight/sumhh
bysort moneypy: egen meankwh=total(weight*kwh)

 collapse meankwh,by(moneypy)


rename moneypy hhfaminc_2005
sort hhfaminc_2005 
quietly destring hhfaminc_2005,replace

quietly merge 1:1 hhfaminc_2005 using  `i'_cps_2005

quietly drop if meankwh==.
quietly drop if hh_count==.

*sum kwh // this sum is not estimated/ estimated is by weights, this is just summary
quietly sum meankwh[aweight=weight]
quietly return list
dis `i',r(mean) 

}


***standard errors

*****====standard errors====*****

clear
use recs_2005

sort moneypy division

egen hhgroup=group(moneypy)

reg kwh i.hhgroup [weight=nweight], nocons 
estimates store results1
esttab results1, b(3) se(3) scalars(r2) star( * 0.10 ** 0.05 *** 0.01) 


matrix B=e(V)

local dim (`= rowsof(A)',`=colsof(A)') 

di "`dim'" 

forval i=1990/2020 {
quietly use `i'_cps_2005,clear
quietly mkmat weight,matrix(A)
quietly matrix SD=A'*B*A
quietly matrix list SD
quietly scalar z=SD[1,1]
di sqrt(z)
}


**** predicted growth


**** generate cps_division count

clear
use income_b2005.dta
bysort year: egen sumhh=total(asecfwt)
bysort year: gen weight=asecfwt/sumhh
bysort year: egen meanincome=total(weight*hhincome_2005)

collapse meanincome,by(year)

*drop if year<2005 //previous after 2005 now all

gen grate=meanincome/meanincome[_n-1]
replace grate=1 if missing(grate)



gen lngrate = ln(grate)
gen lntotal = sum(lngrate)
gen double grateproduct = exp(lntotal)
drop grate
rename grateproduct grate
drop lngrate
drop lntotal
drop meanincome

save grate_05b_1.dta,replace


clear
use income_b2005.dta
merge m:1 year using grate_05b_1 
drop _merge
sort year

gen nn_income=hhincome_2005*grate
order hhincome nn_income grate

save income_b2005_grate.dta,replace
** 1990 -2020
clear

forvalues num=1990/2020 {
use income_b2005_grate.dta
	*** new income groups

gen ngroup = 1 if nn_income<2500
	replace ngroup = 2 if nn_income>=2500 & nn_income<5000
	replace ngroup = 3 if nn_income>=5000 & nn_income<7500
	
	replace ngroup = 4 if nn_income>=7500 & nn_income<10000
	replace ngroup = 5 if nn_income>=10000 & nn_income<15000
	replace ngroup = 6 if nn_income>=15000 & nn_income<20000
	replace ngroup = 7 if nn_income>=20000 & nn_income<25000
	replace ngroup = 8 if nn_income>=25000 & nn_income<30000
	
	replace ngroup = 9 if nn_income>=30000 & nn_income<35000
	replace ngroup = 10 if nn_income>=35000 & nn_income<40000
	
	replace ngroup = 11 if nn_income>=40000 & nn_income<45000
	replace ngroup = 12 if nn_income>=45000 & nn_income<50000
	replace ngroup = 13 if nn_income>=50000 & nn_income<55000
	replace ngroup = 14 if nn_income>=55000 & nn_income<60000
	replace ngroup = 15 if nn_income>=60000 & nn_income<65000
	replace ngroup = 16 if nn_income>=65000 & nn_income<70000
	replace ngroup = 17 if nn_income>=70000 & nn_income<75000
	replace ngroup = 18 if nn_income>=75000 & nn_income<80000
	replace ngroup = 19 if nn_income>=80000 & nn_income<85000
	replace ngroup = 20 if nn_income>=85000 & nn_income<90000
	replace ngroup = 21 if nn_income>=90000 & nn_income<95000
	
	replace ngroup = 22 if nn_income>=95000 & nn_income<100000
	replace ngroup = 23 if nn_income>=100000 & nn_income<120000
	replace ngroup = 24 if nn_income>=120000


keep if year==`num'

collapse (sum) asecwth, by(ngroup)
rename  asecwth hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save grate`num'_cps_05,replace
}


*** predictefd kwh

clear
forval i=1990/2020  {
clear
use recs_2005 // use 1993 year only
** get weighted kwh
bysort moneypy: egen sumhh=total(nweight)
bysort moneypy: gen weight=nweight/sumhh
bysort moneypy: egen meankwh=total(weight*kwh)

collapse meankwh,by(moneypy)

rename moneypy ngroup
quietly sort ngroup 
quietly destring ngroup,replace

quietly merge 1:1 ngroup using  grate`i'_cps_05

quietly drop if meankwh==.
quietly drop if hh_count==.

*sum kwh // this sum is not estimated/ estimated is by weights, this is just summary
quietly sum meankwh[aweight=weight]
quietly return list
dis `i',r(mean)

}


*** standard errors of predicted growth
clear
use recs_2005

sort moneypy division

egen hhgroup=group(moneypy )

reg kwh i.hhgroup [weight=nweight], nocons 
estimates store results1
esttab results1, b(3) se(3) scalars(r2) star( * 0.10 ** 0.05 *** 0.01) 


matrix B=e(V)

local dim (`= rowsof(A)',`=colsof(A)') 

di "`dim'" 

local dim (`= rowsof(B)',`=colsof(B)') 

di "`dim'" 

forval i=1990/2020 {
quietly use grate`i'_cps_05,clear
quietly mkmat weight,matrix(A)
quietly matrix SD=A'*B*A
quietly matrix list SD
quietly scalar z=SD[1,1]
di sqrt(z)
}
